//  ******************************************************************
//  GKC
//
//  Created by Martell on 2012-01-09.
//  Copyright (c) 2012. All rights reserved.
//  Comments:  This is a near copy of Sideeks Golden King Crab Code
//             that I developed to get a better understanding of 
//             his model.
//  TODO: add parameter ivals, bounds and phz to control file.
//  ******************************************************************


DATA_SECTION
	
	init_adstring data_file;
	init_adstring control_file;
	
	// read in the data file
	!! ad_comm::change_datafile_name(data_file);
	
	// model dimensions
	init_int commstyr;  //commercial fishery start year of input catch at length
	init_int commendyr; //commercial fishery & observer data end year of input catch or CPUE at length
	init_int nbins;     // number of length bins for all types of data
	//
	init_matrix retdcpue(commstyr,commendyr,1,nbins);       // total observer retained cpue by length and year
	init_matrix discdcpue(commstyr+5,commendyr,1,nbins);    // total observer discard cpue by length and year
	init_matrix retdcatch(commstyr,commendyr,1,nbins);      // total commercial retained catch by length and year
	init_matrix discdcatch(commstyr+5,commendyr,1,nbins);   // total commercial discarded catch by length and year, 20% mortality aplied
	//
	init_ivector length_bins(1,nbins);                      // mid length of each bin in mm CL
	init_vector  retweight(1,nbins);                        // weight at length in gram
	init_vector  eltime1(commstyr,commendyr);               // elapsed time from survey date to mid fishery time
	init_vector  retcatcheffsample(commstyr,commendyr);     // Retained catch composition effective sample size for 1985 to 2010
	init_vector  discdcatcheffsample(commstyr,commendyr);   // Discard catch composition effective sample size for 1985 to 2010 
	init_vector  gdiscdcatcheffsample(commstyr,commendyr);  // Ground fish Discard catch composition effective sample size for 1985 to 2010
	init_vector  retdcpuevar(commstyr,commendyr);           // Observer retained cpue variance 1985 to 2010
	init_vector  discdcpuevar(commstyr+5,commendyr);        // Observer discard cpue variance 1990 to 2010
	init_matrix  gdiscdcatch(commstyr+10,commendyr,1,nbins);// groundfish discarded catch by length and year
	//
	init_int eof;
	!! if(eof != 999){cout<<"Error reading data, aborting run."<< endl;exit(1);}
	
	
	
	// read in the control file elements
	!! ad_comm::change_datafile_name(control_file);
	init_number alphar;                                     // male mean length of recruits
	init_int recbin;                                        // number of length bins for recruit distribution
	init_int sel_phase;                                     // selectivity phase
	init_vector like_wght(1,18);                            // likelihood weights
	init_number m_disc;                                     // pot fishery handling mortality
	init_number m_gdisc;                                    // groundfish discard
	//
	init_int eocf;
	
	!! if(eocf != 999){cout<<"Error reading control file, aborting run."<<endl; exit(1);}
	!! cout<<"End of reading control file"<<endl;
	!! cout<<alphar<<endl;
	!! cout<<recbin<<endl;
	!! cout<<sel_phase<<endl;
	!! cout<<like_wght<<endl;
	!! cout<<m_disc<<endl;
	!! cout<<m_gdisc<<endl;
	
	//
	vector obsretdcatchN(commstyr,commendyr);				// annual observed retained catch number
	vector obsretdcatchB(commstyr,commendyr);				// annual observed retained catch biomass
	vector obstotalcatchN(commstyr+5,commendyr);
	vector obstotalcatchB(commstyr+5,commendyr);
	vector obsdiscdcatchN(commstyr+5,commendyr); 			// annual observed discarded catch number
	vector obsdiscdcatchB(commstyr+5,commendyr); 			// annual observed discarded catch biomass
	vector obsgdiscdcatchN(commstyr+10,commendyr);  		// annual observed groundfish bycatch number
	vector obsgdiscdcatchB(commstyr+10,commendyr);      	// annual observed groundfish bycatch biomass
	vector obstretdcpue(commstyr,commendyr);
	vector obstdiscdcpue(commstyr+5,commendyr);
	matrix obsretdcpue(commstyr,commendyr,1,nbins);
	matrix obsdiscdcpue(commstyr+5,commendyr,1,nbins);
	vector obstotalcpue(commstyr+5,commendyr);
	matrix obstotcpue(commstyr+5,commendyr,1,nbins);
	matrix obsretlencomp(commstyr,commendyr,1,nbins);
	matrix obstotallencomp(commstyr+5,commendyr,1,nbins);
	matrix obsdiscdlencomp(commstyr+5,commendyr,1,nbins);
	matrix obsgdiscdlencomp(commstyr+10,commendyr,1,nbins); 
	//
	
	LOC_CALCS
		// calculations only need to be done once, here you can allocate fixed parameter values
		int t,l,j;
		// changed the effective sample size  SHOULD GO IN DATA FILE
		// 
		retcatcheffsample    = 250.;    //50
		discdcatcheffsample  = 10.;     //50
		gdiscdcatcheffsample = 10.;     //just before used 200, original10
		//
		retdcatch   /= 1000000.;    				// retained catch in millions of crabs
		discdcatch   = discdcatch*m_disc/200000.;   // discard catch in millions of crabs, 
													// hm 20% already applied but allowed for any hm changes
		gdiscdcatch /= 1000000.;   					// groundfish discard catch in millions of crabs
		//
		// Observed retained,pot and groundfish discard, and total pot catch by year in millions of crabs/metric tons
		//
		obsretdcatchN.initialize();
		obsdiscdcatchN.initialize();
		obsretdcatchB.initialize();
		obsdiscdcatchB.initialize();
		obsgdiscdcatchN.initialize();
		obsgdiscdcatchB.initialize();
		obstotalcatchN.initialize();
		obstotalcatchB.initialize();
		//
		for(t = commstyr; t <= commendyr; t++)
		{
			for(l = 1; l <= nbins; l++)
			{
				obsretdcatchN(t) += retdcatch(t,l);
				obsretdcatchB(t) += retdcatch(t,l)*retweight(l);                        // in metric tons
				if (t>(commstyr+4))obsdiscdcatchN(t)  += discdcatch(t,l);
				if (t>(commstyr+4))obsdiscdcatchB(t)  += discdcatch(t,l)*retweight(l);  // in metric tons
				if (t>(commstyr+9))obsgdiscdcatchN(t) += gdiscdcatch(t,l);
				if (t>(commstyr+9))obsgdiscdcatchB(t) += gdiscdcatch(t,l)*retweight(l); // in metric tons
			}
			if (t>(commstyr+4)) obstotalcatchN(t) = obsretdcatchN(t)+obsdiscdcatchN(t);
			if (t>(commstyr+4)) obstotalcatchB(t) = obsretdcatchB(t)+obsdiscdcatchB(t);
		}
		//
		obsretlencomp.initialize();
		obsdiscdlencomp.initialize();
		obstotallencomp.initialize();
		obsgdiscdlencomp.initialize();
		//
		for(t = commstyr; t <= commendyr; t++)
		{
			for(l = 1; l <= nbins; l++)
			{
				obsretlencomp(t,l) = retdcatch(t,l)/obsretdcatchN(t);
				if (t>(commstyr+4))obsdiscdlencomp(t,l)  = discdcatch(t,l)/obsdiscdcatchN(t);
				if (t>(commstyr+4))obstotallencomp(t,l)  = (retdcatch(t,l)+discdcatch(t,l))/obstotalcatchN(t);
				if (t>(commstyr+9))obsgdiscdlencomp(t,l) = gdiscdcatch(t,l)/obsgdiscdcatchN(t);
			}
		}
		//
		// Observed retained, discard, and total CPUE by year in number of crabs
		// The values in 1993 and 1994 are set to zero for E174W because no observer sampling
		//
		obsretdcpue.initialize();
		obstretdcpue.initialize();
		obsdiscdcpue.initialize();
		obstdiscdcpue.initialize();
		obstotcpue.initialize();
		obstotalcpue.initialize();
		//
		for(t = commstyr; t <= commendyr; t++)
		{
			//cout<<"Martell was here"<<t<<endl;
			for(l = 1; l <= nbins; l++)
			{
				obsretdcpue(t,l) = retdcpue(t,l);
				obstretdcpue(t) += retdcpue(t,l);
				if (t>(commstyr+4))obsdiscdcpue(t,l) = discdcpue(t,l);
				if (t>(commstyr+4))obstdiscdcpue(t) += discdcpue(t,l);
				if (t>(commstyr+4))obstotcpue(t,l)   = retdcpue(t,l)+discdcpue(t,l);
			}
			if (t>(commstyr+4))obstotalcpue(t)=obstretdcpue(t)+obstdiscdcpue(t);
		}
	END_CALCS

INITIALIZATION_SECTION
	// initial values for estimating parameters
	// this should all go into a control file, or the control file so its not hard wired.
	log_mean_rec 0.1    
	log_mean_Fpot -2.3  
	log_mean_Fground -4.6  
	log_newsh1 1.61   
	log_betar -1.273     
	log_beta  -0.548     
	a 0.05               
	log_b 4.935           
	c2 0.02
	log_d2 4.8
	c3 0.02
	log_d3 4.8
	aa2 0.02
	log_bb2  4.89
	aa3 0.02
	log_bb3   4.89             
	aa4 0.02
	log_bb4   4.89
	logq2 -13.
	logq3 -13.
	ret_var 0.05
	discd_var 0.05


	
PARAMETER_SECTION
	//this section's variables are global, can use C++ codes within LOCAL_CALCS  ending with END_CALCS
	//SJDM, this is a bit scary with all the hard wired bounds. Should all go into a control file.
	//Idea would be to have one source code file and data files and control files for each area.
	
	
	init_bounded_dev_vector rec_dev(commstyr+1,commendyr+1,-15.,15.,3);  
	init_bounded_dev_vector Fpot_dev(commstyr,commendyr,-15.,15.,2);    
	init_bounded_dev_vector Fground_dev(commstyr+10,commendyr,-15.,15.,2); 
	init_bounded_number     a(0.01,0.75,2);                      // molt probability slope
	init_bounded_number     log_b(3.98,5.5,2);                   // molt probability L50
	init_bounded_number     log_beta(-0.7,2.0,2);                // growth increment earlier logbeta(-3.6,2.,1)
	init_bounded_number     c2(0.01,0.8,sel_phase);              // fishery total selectivity slope for <=2004
	init_bounded_number     log_d2(3.98,4.98,sel_phase);         // fishery total selectivity L50
	init_bounded_number     c3(0.01,0.8,sel_phase);              // fishery total selectivity slope for 2005 - 
	init_bounded_number     log_d3(3.98,5.5,sel_phase);          // fishery total selectivity L50
	init_bounded_number     aa2(0.01,1.5,sel_phase);             // retained selectivity slope for <=2004
	init_bounded_number     log_bb2(4.85,4.98,sel_phase);        // retained selectivity L50
	init_bounded_number     aa3(0.01,3.0,sel_phase);             // retained selectivity slope for 2005 - 
	init_bounded_number     log_bb3(4.25,4.98,sel_phase);        // retained selectivity L50
	init_bounded_number     aa4(0.01,1.5,sel_phase);             
	init_bounded_number     log_bb4(3.98,5.5,sel_phase);         // groundfish bycatch selectivity L50
	// init_bounded_number log_betar(-4.6,1.2,3);                // recruitment betar
	init_bounded_number     log_betar(-4.6,3.2,3);               // recruitment betar
	init_bounded_number     logq2(-25.,-0.01,sel_phase+1);       // pot fishery catchability in log form for <=2004
	init_bounded_number     logq3(-25.,-0.01,sel_phase+1);       // pot fishery catchability in log form for 2005 -
	init_bounded_number     log_newsh1(0.01,10.,1);              // start year newshell abundance
	// init_bounded_number log_oldsh1(0.01,6.,1);                // start year oldshell abundance
	init_bounded_number     log_mean_rec(0.01,3.,2);             // log mean recruit 
	init_bounded_number     log_mean_Fpot(-15,-0.35,1);          // log mean F in pot fishery
	init_bounded_number     log_mean_Fground(-15.,-1.6,1);       // log mean F in groundfish fishery bycatch
	init_bounded_dev_vector qq_dev(1,2,-20.,20.,sel_phase+2);    // legal discard
	init_number             mean_qq(sel_phase);                  // legal discard  
	init_number             M(-1);                               // natural mortality
	!! M=0.18;
	init_bounded_vector     alphaN(1,nbins,-5.,5.,1);            // new shell initial size diatribution 
	init_bounded_number     ret_var(0.01,1.25,sel_phase+1);      // ret cpue extra variance
	init_bounded_number     discd_var(0.01,1.25,sel_phase+1);    // pot discard cpue extra variance
	init_bounded_number     Ftemp(0.,0.75);                      // Fofl iteration
	
	// special variables
	sdreport_vector     mat_biomass(commstyr,commendyr+1);
	sdreport_vector     legal_biomass(commstyr,commendyr+1);
	likeprof_number MMB_lprof;
	// likeprof_number predtotofl_lprof;
	// likeprof_number predretofl_lprof;
	// likeprof_number M_lprof;
	//
	//
	objective_function_value f;
	//
	//
	number surv1;
	number surv2;
	number surv3;
	number surv4;
	number surv5;
	number preddiscdcatchBB;
	number predgdiscdcatchBB;    //
	number final_F;
	number retain_harvest;
	number total_harvest;
	number meanmatbiomass;
	number like_retdcatchB;
	number like_discdcatchB;
	number like_gdiscdcatchB;    //
	number like_rec_dev;
	number like_retlencomp;
	number like_discdlencomp;
	number like_gdiscdlencomp;    //
	number like_retcpue;
	number like_discdcpue;
	number like_F;
	number like_Mpenalty;
	number like_initialLF;
	number like_smoothLFQ;
	number like_gF;
	number like_qq;
	number like_molt50;
	number like_meanF;
	number like_beta;
	number like_meanq;
	number like_finalF;
	number qtemp;
	number recsum;
	number alpha;
	number alpha_rec;
	number temp5;
	number temp6;
	number temp7;
	number temp8;
	number inffmort;
	number ffmort;
	number maxsel_totsel;
	number maxsel_retsel;
	number maxsel_gsel;
	number grow_incr;
	number mean_gfmort;
	number predofl;
	number t5;
	number betal;


	vector totalttemp(commstyr,commendyr);
	vector totalttempB(commstyr,commendyr);
	vector retcatchsum(commstyr,commendyr);
	vector discdcatchsum(commstyr,commendyr);          //
	vector gdiscdcatchsum(commstyr,commendyr);       //
	vector predttotcatchN(commstyr,commendyr);
	vector predretdcatchB(commstyr,commendyr);
	vector preddiscdcatchB(commstyr,commendyr);
	vector predtgdiscdcatchN(commstyr,commendyr);         //
	vector predgdiscdcatchB(commstyr,commendyr);          //
	vector predttotalcpue(commstyr,commendyr);
	vector predtretcpue(commstyr,commendyr);
	vector predtdiscdcpue(commstyr,commendyr);
	vector expalphaN(1,nbins);
	vector first_propN(1,nbins);
	vector first_propO(1,nbins);
	vector legalabund(commstyr,commendyr+1);
	vector matabund(commstyr,commendyr+1);
	//vector legal_biomass(commstyr,commendyr+1);
	//vector mat_biomass(commstyr,commendyr+1);
	vector harvestrate(commstyr,commendyr);
	vector rec_len(1,nbins);
	vector qq(commstyr,commendyr+1);
	vector qqtemp(1,2);
	vector gselect(1,nbins);
	vector molt(1,nbins);
	vector sum_len(1,nbins);
	vector devia(1,51);
	vector RR(commstyr,commendyr+1);
	vector gfmort_tem(commstyr+10,commendyr);
	vector gfmort_final(commstyr,commendyr);
	vector t11(1,60);
	vector tt1(1,nbins);
	vector tt2(1,nbins);
	vector explegalB(commstyr,commendyr+1);
	vector legalcatch(commstyr,commendyr);
	
	
	matrix newsh(commstyr,commendyr+1,1,nbins);
	matrix oldsh(commstyr,commendyr+1,1,nbins);
	matrix discdtemp(commstyr,commendyr,1,nbins);
	matrix gdiscdtemp(commstyr,commendyr,1,nbins);    //
	matrix retcatchlen(commstyr,commendyr,1,nbins);
	matrix retsigmasquare(commstyr,commendyr,1,nbins);
	matrix discdsigmasquare(commstyr+5,commendyr,1,nbins);    //
	matrix gdiscdsigmasquare(commstyr+10,commendyr,1,nbins);  //
	matrix discdcatchlen(commstyr,commendyr,1,nbins);    //
	matrix gdiscdcatchlen(commstyr,commendyr,1,nbins);   //
	matrix predretcatchlen(commstyr,commendyr,1,nbins);
	matrix preddiscdcatchlen(commstyr,commendyr,1,nbins);
	matrix predgdiscdcatchlen(commstyr,commendyr,1,nbins);  //
	matrix predtotalcpue(commstyr,commendyr,1,nbins);
	matrix preddiscdcpue(commstyr,commendyr,1,nbins);
	matrix predretcpue(commstyr,commendyr,1,nbins);
	matrix predtotcatchN(commstyr,commendyr+1,1,nbins);
	matrix predretdcatchN(commstyr,commendyr+1,1,nbins);
	matrix preddiscdcatchN(commstyr,commendyr+1,1,nbins);
	matrix predgdiscdcatchN(commstyr,commendyr+1,1,nbins);  //
	matrix len_len(1,nbins,1,nbins);
	matrix totalselect(commstyr,commendyr+1,1,nbins);
	matrix retselect(commstyr,commendyr+1,1,nbins);
	matrix ret_stdresid(commstyr,commendyr,1,nbins);
	matrix discd_stdresid(commstyr+5,commendyr,1,nbins);
	matrix gdiscd_stdresid(commstyr+10,commendyr,1,nbins);
	matrix tt3(1,nbins,1,60);
	matrix ttt(1,nbins,1,60);
	
	
	
PROCEDURE_SECTION
	
	/*________ MAIN ROUTINES ________*/
	recruit();               // recruit ditribution  (rec_len)
	molting();               // molt probabilty
	grow();                  // this has growth matrix (len_len)
	totalselectivity();      // total selectivity
	retselectivity();        // retained selectivity
	gselectivity();          // groundfish discard selectivity
	gmortality();            // groundfish bycatch fishing mortality
	population_model();	     // initialize and run the population model, retained and discard catch
	calc_predicted_cpue();   // calculate the predicted cpue
	wtf();				     // this needs some major work, a lot of code here that needs to be parsed.
	
	objective_function();  //ALSO NEEDS WORK!
	/*_______________________________*/

FUNCTION population_model
	int t,i,l,j,ll;
	qqtemp.initialize();
	for(j = 1; j <=2; j++)
	{
		qqtemp(j)=mfexp(mean_qq+qq_dev(j)*0.25)/(1.+mfexp(mean_qq+qq_dev(j)*0.25));
	} 
	//
	surv1 = mfexp(-M);
	//
	newsh.initialize();
	oldsh.initialize();
	predtotcatchN.initialize();
	predttotcatchN.initialize();
	predretdcatchN.initialize();
	predretdcatchB.initialize();
	preddiscdcatchN.initialize();
	preddiscdcatchB.initialize();
	predgdiscdcatchN.initialize();
	predgdiscdcatchB.initialize();
	predtgdiscdcatchN.initialize();
	legalcatch.initialize();
	//
	expalphaN.initialize();
	//   expalphaO.initialize();
	RR.initialize();
	ffmort.initialize();
	//
	// Initial newshell and oldshell relative size frequency proportion parameters
	//
	for(l = 1; l <= nbins; l++)
	{
		expalphaN(l) = mfexp(alphaN(l));
	}
	first_propN.initialize();
	first_propO.initialize();
	for(l = 1; l <= nbins; l++)
	{
		if(sum(expalphaN)>0.)
			first_propN(l) = expalphaN(l)/sum(expalphaN);
		if(first_propN(l)<=0.)
			first_propN(l) = 0.0001;
		first_propO(l)  = 0.;
	}
	//
	// Population dynamics
	//
	for(t = commstyr; t <= commendyr; t++)
	{
		qq(t)=1.;
		if(t==1991)
			qq(t)=qqtemp(1);
		if(t==1994)
			qq(t)=qqtemp(2);
		
		surv2   = mfexp((eltime1(t)-1.)*M);
		surv3   = mfexp(-eltime1(t)*M);
		ffmort  = mfexp(log_mean_Fpot+Fpot_dev(t));
		RR(t+1) = mfexp(log_mean_rec+rec_dev(t+1));
		if(t==commendyr)
			RR(commendyr+1)=(RR(commendyr-2)+ RR(commendyr-1)+RR(commendyr))/3.;
		//
		if(t == commstyr)
		{
			for(l = 1; l <= nbins; l++)
			{
				newsh(t,l) = mfexp(log_newsh1)*first_propN(l); // in millions, proportion parameter is estimated in the model
				oldsh(t,l) = 0.;
			}
		}
		//


		// Predicted retained and discard catch
		// SJDM: I don't think the series of if statements in the following loop should be there b/c the function
		// is disconitnous other wise.  Currently testing to see if any are true, and forcing exit.
		// I checked each one, with the adanom8510aps data, and did not have any exits.  I've left in the exits
		// for checking with other data files.
		for(l = 1; l <= nbins; l++)
		{
			predtotcatchN(t,l)    = (newsh(t,l) + oldsh(t,l)) * surv3*(1.-mfexp(-ffmort*totalselect(t,l)-gfmort_final(t)*gselect(l)))
									* ffmort*totalselect(t,l)/(ffmort*totalselect(t,l)+gfmort_final(t)*gselect(l));  // millions of crabs
			if(predtotcatchN(t,l)<=0.){
				cout<<predtotcatchN(t,l)<<endl;
				exit(1);
				//if this is true, then the function is not continuous!
				predtotcatchN(t,l)=0.0;
			}
			 
			predttotcatchN(t)    += predtotcatchN(t,l);
			if(predttotcatchN(t)<=0.){
				cout<<predttotcatchN(t,l)<<endl;
				exit(1);
				predttotcatchN(t)=0.0000001;
			}
			
			predretdcatchN(t,l)   = (newsh(t,l) + oldsh(t,l)) * surv3*(1.-mfexp(-ffmort*totalselect(t,l)-gfmort_final(t)*gselect(l)))
									*ffmort*totalselect(t,l)*qq(t)*retselect(t,l)/(ffmort*totalselect(t,l)+gfmort_final(t)*gselect(l));
			if(predretdcatchN(t,l)<=0.){
				cout<<predretdcatchN(t,1)<<endl;
				exit(1);
				predretdcatchN(t,l)=0.0;
			}
			
			predretdcatchB(t) += predretdcatchN(t,l)*retweight(l); // metric tons
			if(predretdcatchB(t)<=0.){
				cout<<predretdcatchB(t)<<endl;
				exit(1);
				predretdcatchB(t)=0.00000001;
			}
			
			preddiscdcatchN(t,l)  = m_disc*(predtotcatchN(t,l)-predretdcatchN(t,l));  // handling mortality applied
			//cout<<preddiscdcatchN(t,l)<<endl;
			if(preddiscdcatchN(t,l)<=0.){
				//cout<<m_disc<<endl; //the value of m_disc is 0.2 (20% handling mortality)
				if(preddiscdcatchN(t,l)<0.){ 
					/*SJDM This one is always 0 when predtotcatchN-predretdcatchN =0 */
					cout<<preddiscdcatchN(t,l)<<endl;
					exit(1);
				}
				preddiscdcatchN(t,l)=0.0;
			}
			
			preddiscdcatchB(t) += preddiscdcatchN(t,l)*retweight(l);
			if(preddiscdcatchB(t)<=0.){
				cout<<preddiscdcatchB(t)<<endl;
				exit(1);
				preddiscdcatchB(t)=0.00000001;
			}
			
			predgdiscdcatchN(t,l) = m_gdisc*(newsh(t,l) + oldsh(t,l)) * surv3*(1.-mfexp(-ffmort*totalselect(t,l)-gfmort_final(t)*gselect(l)))
									*gfmort_final(t)*gselect(l)/(ffmort*totalselect(t,l)+gfmort_final(t)*gselect(l));  // millions of crabs
			if(predgdiscdcatchN(t,l)<=0.){
				cout<<predgdiscdcatchN(t,l)<<endl;
				exit(1);
				predgdiscdcatchN(t,l)=0.0;
			}
			
			predtgdiscdcatchN(t) += predgdiscdcatchN(t,l);
			if(predtgdiscdcatchN(t)<=0.){
				cout<<predtgdiscdcatchN(t)<<endl;
				exit(1);
				predtgdiscdcatchN(t)=0.0000001;
			}
			
			predgdiscdcatchB(t) += predgdiscdcatchN(t,l)*retweight(l); // metric tons
			if(predgdiscdcatchB(t)<=0.){
				cout<<predgdiscdcatchB(t)<<endl;
				exit(1);
				predgdiscdcatchB(t)=0.0000001;
			}
			//many hard wired switches in this code, it would be nice to have them accessible w/o having to compile
			if(l>7)legalcatch(t) += predretdcatchN(t,l) + preddiscdcatchN(t,l);
		}
		//
		for(l = 1; l <= nbins; l++)
		{
			//
			for(ll = 1; ll <= l; ll++)
			{
				dvariable yyy = ((newsh(t,ll) + oldsh(t,ll))*surv1 - (predretdcatchN(t,ll) + preddiscdcatchN(t,ll)+predgdiscdcatchN(t,ll))*surv2)*molt(ll)*len_len(ll,l);
				newsh(t+1,l) += yyy;                                              // newshell abundance millions
				//   if (current_phase()>7) cout<<l<<" yyy = "<<yyy<<"  newsh = "<<newsh(t+1,l)<<endl;
			}
			newsh(t+1,l) += RR(t+1)*rec_len(l);
			//   if (current_phase()>7) cout<<l<< " rr = "<<RR(t+1)*rec_len(l)<<"  newsh = "<<newsh(t+1,l)<<endl;

			if(newsh(t+1,l)<=0.000001)newsh(t+1,l) = 0.000001;
			dvariable xxx = ((newsh(t,l) + oldsh(t,l))*surv1 - (predretdcatchN(t,l) + preddiscdcatchN(t,l)+predgdiscdcatchN(t,l))*surv2)*(1.0 - molt(l));
			oldsh(t+1,l) = xxx;                                                // oldshell abundance millions
			if(oldsh(t+1,l)<=0.000001)oldsh(t+1,l) = 0.000001;
		}
	}

FUNCTION calc_predicted_cpue
	int t,i,l,j,ll;
	//
	// Predicted total, retained and discard CPUE in the fishery by length and year
	//
	predtotalcpue.initialize();
	predretcpue.initialize();
	preddiscdcpue.initialize();
	for(t = commstyr; t <= commendyr; t++)
	{
		ffmort=mfexp(log_mean_Fpot+Fpot_dev(t));
		qq(t)=1.;
		if(t==1991)qq(t)=qqtemp(1);
		if(t==1994)qq(t)=qqtemp(2);
		surv3=mfexp(-eltime1(t)*M);
		if(t <= 1997)qtemp = logq2;
		if(t >= 1998 && t <= 2004)qtemp = logq2;
		if(t > 2004)qtemp = logq3;
		for(l = 1; l <= nbins; l++)
		{
			predtotalcpue(t,l)  = 1000000.*mfexp(qtemp)*totalselect(t,l)*(newsh(t,l) + oldsh(t,l) - 0.5*(predretdcatchN(t,l)+preddiscdcatchN(t,l)+predgdiscdcatchN(t,l)))*surv3; // number of crabs
			predretcpue(t,l)    = 1000000.*mfexp(qtemp)*totalselect(t,l)*qq(t)*retselect(t,l)*(newsh(t,l) + oldsh(t,l) - 0.5*(predretdcatchN(t,l)+preddiscdcatchN(t,l)+predgdiscdcatchN(t,l)))*surv3;
			preddiscdcpue(t,l)  = predtotalcpue(t,l)-predretcpue(t,l);    // no handling mortality applied
			if(predtotalcpue(t,l)<=0.0){
				exit(1);
				predtotalcpue(t,l) = 0.0001;
			}
			if(predretcpue(t,l)<=0.0){
				exit(1);
				predretcpue(t,l)   = 0.0001;
			}
			if(preddiscdcpue(t,l)<=0.0){
				//cout<<t<<endl;
				//exit(1);  //SJDM: Is this an issue? The code is not differentiable here.
				preddiscdcpue(t,l) = 0.0001;
			}
		}
	}
	
	predttotalcpue.initialize();
	predtretcpue.initialize();
	predtdiscdcpue.initialize();
	for(t = commstyr; t <= commendyr; t++)
	{
		for(l = 1; l <= nbins; l++)
		{
			predttotalcpue(t) +=predtotalcpue(t,l);
			predtretcpue(t) +=predretcpue(t,l);
			predtdiscdcpue(t) += preddiscdcpue(t,l);
		}
	}
	//
	//   Predicted retained and discard catch length composition by year in the fishery
	//
	retcatchsum.initialize();
	discdcatchsum.initialize();
	gdiscdcatchsum.initialize();
	retcatchlen.initialize();
	discdcatchlen.initialize();
	gdiscdcatchlen.initialize();
	//
	for(t = commstyr; t <= commendyr; t++)
	{
		for(l = 1; l <= nbins; l++)
		{
			retcatchlen(t,l) = predretdcatchN(t,l);
			retcatchsum(t)  += predretdcatchN(t,l);
			discdcatchlen(t,l) = preddiscdcatchN(t,l);
			discdcatchsum(t) +=  preddiscdcatchN(t,l);
			gdiscdcatchlen(t,l) = predgdiscdcatchN(t,l);
			gdiscdcatchsum(t) +=  predgdiscdcatchN(t,l);
		}
	}
	predretcatchlen.initialize();
	preddiscdcatchlen.initialize();
	predgdiscdcatchlen.initialize();
	retsigmasquare.initialize();
	discdsigmasquare.initialize();
	gdiscdsigmasquare.initialize();
	for(t = commstyr; t <= commendyr; t++)
		{
		for(l = 1; l <= nbins; l++)
		{
			if(retcatchsum(t)>0.)predretcatchlen(t,l) = retcatchlen(t,l)/retcatchsum(t);
			if(predretcatchlen(t,l)<=0.)predretcatchlen(t,l)=0.0001;
			if(discdcatchsum(t)>0.)preddiscdcatchlen(t,l) = discdcatchlen(t,l)/discdcatchsum(t);
			if(preddiscdcatchlen(t,l)<=0.)preddiscdcatchlen(t,l)=0.0001;
			if(gdiscdcatchsum(t)>0.)predgdiscdcatchlen(t,l) = gdiscdcatchlen(t,l)/gdiscdcatchsum(t);
			if(predgdiscdcatchlen(t,l)<=0.)predgdiscdcatchlen(t,l)=0.0001;
			//
			//   Length composition variance estimation using multinomial formula:
			//
			retsigmasquare(t,l)= ((1.0-obsretlencomp(t,l))*obsretlencomp(t,l)+0.1/double(nbins))/(retcatcheffsample(t));
			if(t>(commstyr+4))discdsigmasquare(t,l)= ((1.0-obsdiscdlencomp(t,l))*obsdiscdlencomp(t,l)+0.1/double(nbins))/(discdcatcheffsample(t));
			if(t>(commstyr+9))gdiscdsigmasquare(t,l)= ((1.0-obsgdiscdlencomp(t,l))*obsgdiscdlencomp(t,l)+0.1/double(nbins))/(gdiscdcatcheffsample(t));
			//
		}
	}
	//
	ret_stdresid.initialize();
	discd_stdresid.initialize();
	gdiscd_stdresid.initialize();
	for(t = commstyr; t <= commendyr; t++)
	{
		for(l = 1; l <= nbins; l++)
		{
			ret_stdresid(t,l)=(obsretlencomp(t,l)-predretcatchlen(t,l))/sqrt(2.*retsigmasquare(t,l));
			if(t>(commstyr+4))discd_stdresid(t,l)=(obsdiscdlencomp(t,l)-preddiscdcatchlen(t,l))/sqrt(2.*discdsigmasquare(t,l));
			if(t>(commstyr+9))gdiscd_stdresid(t,l)=(obsgdiscdlencomp(t,l)-predgdiscdcatchlen(t,l))/sqrt(2.*gdiscdsigmasquare(t,l));
		}
	}
	//
	// Predicted legal male abundance in number (millions) and weight (t) by year at survey time
	//
	legalabund.initialize();
	legal_biomass.initialize();
	explegalB.initialize();
	for(t = commstyr;t <= commendyr+1; t++)
	{
		//      surv3=mfexp(-eltime1(t)*M);
		for(l = 8; l <= nbins; l++)        // bin number 8 starts from 136 mm CL, middle length 138 mm CL
		{
			legalabund(t)    += (newsh(t,l) + oldsh(t,l));
			legal_biomass(t) += (newsh(t,l) + oldsh(t,l))*retweight(l);
			explegalB(t)     += (newsh(t,l) + oldsh(t,l))*retselect(t,l)*totalselect(t,l)*retweight(l);
		}
	}
	//






FUNCTION wtf
	//
	int t,i,l,j,ll;
	if(last_phase())
	{
		//////***************************************************
		//
		//Predicted pot fishery harvest rate based on predicted legal size catch divided by legal abundance
		//
		harvestrate.initialize();
		for(t = commstyr; t <= commendyr; t++)
		{
			if(legalabund(t)>0.)
				harvestrate(t) = legalcatch(t)/legalabund(t);
			else
				harvestrate(t) =0.;
		}
		//
		// Predicted mature male (>=121 mmCL) abundance in number (millions) and weight (t) by year on next 15 February upto 
		// final data avialability year
		//
		matabund.initialize();
		mat_biomass.initialize();
		surv4=mfexp(-0.627397*M);
		for(t = commstyr;t <= commendyr; t++)
		{
			surv5 = mfexp((eltime1(t)-0.627397)*M);
			for(l = 5; l <= nbins; l++)        // bin number 5 starts from 121 mm CL, middle length 123 mm CL
			{
				matabund(t) += (newsh(t,l) + oldsh(t,l))*surv4
								-(predretdcatchN(t,l)+preddiscdcatchN(t,l)+predgdiscdcatchN(t,l))*surv5;
				mat_biomass(t) += ((newsh(t,l) + oldsh(t,l))*surv4
								-(predretdcatchN(t,l)+preddiscdcatchN(t,l)+predgdiscdcatchN(t,l))*surv5)*retweight(l);
			}
			if(matabund(t)<=0.)matabund(t)=0.0001;
			if(mat_biomass(t)<=0.)mat_biomass(t)=0.0001;
		}
		meanmatbiomass=sum(mat_biomass)/(commendyr-(commstyr)+1);
		//
		MMB_lprof=meanmatbiomass;
		//
		// Predicted mature male (>=121 mmCL) abundance in number (millions) and weight (t) on next 15 February for the one step forward projection
		// 
		matabund(commendyr+1)=0.;
		mat_biomass(commendyr+1)=0.;
		qq(commendyr+1)=1.;
		//
		for(l = 1; l <= nbins; l++)
		{
			predtotcatchN(commendyr+1,l) = (newsh(commendyr+1,l) + oldsh(commendyr+1,l))*surv3*(1.-mfexp(-Ftemp*totalselect(commendyr+1,l)-gfmort_final(commendyr)*gselect(l)))*Ftemp*totalselect(commendyr+1,l)/(Ftemp*totalselect(commendyr,l)+gfmort_final(commendyr)*gselect(l));  // millions of crabs
			if(predtotcatchN(commendyr+1,l)<=0.)predtotcatchN(commendyr+1,l)=0.0;
			predretdcatchN(commendyr+1,l) = (newsh(commendyr+1,l) + oldsh(commendyr+1,l))*surv3*(1.-mfexp(-Ftemp*totalselect(commendyr+1,l)-gfmort_final(commendyr)*gselect(l)))*Ftemp*totalselect(commendyr+1,l)*qq(commendyr+1)*retselect(commendyr+1,l)/(Ftemp*totalselect(commendyr+1,l)+gfmort_final(commendyr)*gselect(l));
			if(predretdcatchN(commendyr+1,l)<=0.)predretdcatchN(commendyr+1,l)=0.0;
			preddiscdcatchN(commendyr+1,l) = m_disc*(predtotcatchN(commendyr+1,l)-predretdcatchN(commendyr+1,l));  // handling mortality applied
			if(preddiscdcatchN(commendyr+1,l)<=0.)preddiscdcatchN(commendyr+1,l)=0.0;
			predgdiscdcatchN(commendyr+1,l) = m_gdisc*(newsh(commendyr+1,l) + oldsh(commendyr+1,l))*surv3*(1.-mfexp(-Ftemp*totalselect(commendyr+1,l)-gfmort_final(commendyr)*gselect(l)))*gfmort_final(commendyr)*gselect(l)/(Ftemp*totalselect(commendyr+1,l)+gfmort_final(commendyr)*gselect(l));  // millions of crabs
			if(predgdiscdcatchN(commendyr+1,l)<=0.)predgdiscdcatchN(commendyr+1,l)=0.0;
		}
		//
		//SJDM: Question what is this constant 0.627397 for?
		surv5 = mfexp((eltime1(commendyr)-0.627397)*M);
		for(l = 5; l <= nbins; l++)        // bin number 5 starts from 121 mm CL, middle length 123 mm CL
		{
			matabund(commendyr+1) += (newsh(commendyr+1,l) + oldsh(commendyr+1,l))*surv4-(predretdcatchN(commendyr+1,l)+preddiscdcatchN(commendyr+1,l)+predgdiscdcatchN(commendyr+1,l))*surv5;
			mat_biomass(commendyr+1) += ((newsh(commendyr+1,l) + oldsh(commendyr+1,l))*surv4-(predretdcatchN(commendyr+1,l)+preddiscdcatchN(commendyr+1,l)+predgdiscdcatchN(commendyr+1,l))*surv5)*retweight(l);
		}
		if(matabund(commendyr+1)<=0.)matabund(commendyr+1)=0.0001;
		if(mat_biomass(commendyr+1)<=0.)mat_biomass(commendyr+1)=0.0001;
		//
		// predicted F_OFL 
		//
		if(mat_biomass(commendyr+1)>=meanmatbiomass)final_F = M;
		if(mat_biomass(commendyr+1)<meanmatbiomass && mat_biomass(commendyr)>0.25*meanmatbiomass)
		{final_F = M*(mat_biomass(commendyr+1)/meanmatbiomass-0.1)/0.9;}
		if(mat_biomass(commendyr+1)<= 0.25*meanmatbiomass)final_F=0.;
		//
		//  Use the F to get the predicted harvest for 2011/12:
		//
		retain_harvest=0.;
		total_harvest=0.;
		preddiscdcatchBB=0.;
		predgdiscdcatchBB=0.;
		predofl=0.;
		qq(commendyr+1)=1.;
		surv3=mfexp(-eltime1(commendyr)*M);
		for(l = 1; l <= nbins; l++)
		{
			total_harvest += retweight(l)*(newsh(commendyr+1,l) + oldsh(commendyr+1,l))*surv3*(1.-mfexp(-final_F*totalselect(commendyr+1,l)-gfmort_final(commendyr)*gselect(l)))*final_F*totalselect(commendyr+1,l)/(final_F *totalselect(commendyr+1,l)+gfmort_final(commendyr)*gselect(l));
			retain_harvest += retweight(l)*(newsh(commendyr+1,l) + oldsh(commendyr+1,l))*surv3*(1.-mfexp(-final_F*totalselect(commendyr+1,l)-gfmort_final(commendyr)*gselect(l)))*final_F*totalselect(commendyr+1,l)*qq(commendyr+1)*retselect(commendyr+1,l)/(final_F *totalselect(commendyr+1,l)+gfmort_final(commendyr)*gselect(l));
			predgdiscdcatchBB += retweight(l)*m_gdisc*(newsh(commendyr+1,l) + oldsh(commendyr+1,l))*surv3*(1.-mfexp(-final_F*totalselect(commendyr+1,l)-gfmort_final(commendyr)*gselect(l)))*gfmort_final(commendyr)*gselect(l)/(final_F*totalselect(commendyr+1,l)+gfmort_final(commendyr)*gselect(l));
		}
		preddiscdcatchBB = m_disc*(total_harvest-retain_harvest);

		predofl=retain_harvest+predgdiscdcatchBB+preddiscdcatchBB;
		//
		//   predtotofl_lprof=predofl; 
		//   predretofl_lprof=retain_harvest;    
		//
		////******************************************************************
		//
	}                          // last phase end brace




FUNCTION recruit
	////
	int ilen, il2, j, ii;
	dvariable t2, t3, t4;
	rec_len.initialize();
	t11.initialize();
	/*
	The gamma function is truncated and numerically simplified.
	TODO: This could be replaced with my dgamma function, not sure if this is continuous and differentiable.
	*/
	alpha_rec = mfexp(log_betar);
	t4        = alphar/alpha_rec;    // alphar is the mean recruit length
	ii        = (recbin)*5;
	ilen      = 5;
	t2        = 0.0;
	t3        = (t4-1.0)*log(alphar)-t4-t4*log(alpha_rec);
	for (j = 1; j<=ii;j++)
	{
		t5      = ((j-1)+100.0)*1.0;
		t11(j)  = (t4-1.0)*log(t5)-t5/alpha_rec-t4*log(alpha_rec)-t3;
		if (t11(j) < -1.0e30) 	// SJDM, this is not continuous.
			t11(j) = -1.0e30;	
		t11(j)  = mfexp(t11(j));
		t2      = t2 + t11(j);
	}
	for (j=1;j<=recbin-0;j++)
	{
		for (il2=1;il2<=ilen;il2++) 
			rec_len(j) += t11((j-1)*ilen+il2)/t2;
	}

FUNCTION molting
	//assuming a logistic function
	int j;
	for(j = 1; j <= nbins; j++)
	{
		molt(j)=1./(1.+mfexp(a*(double(length_bins(j))-mfexp(log_b))));
	}

FUNCTION grow
	// linear growth distributions using gamma function
	// SJDM, this could be cleaned up with dgamma
	//
	int i, j, ii;
	dvariable tt4, t2;
	len_len.initialize();
	tt3.initialize();
	ttt.initialize();
	betal = mfexp(log_beta);
	for (j=1;j<=nbins;j++)
	{
		t2     = (18.173-0.029*(j*5.0+98.0));
		tt1(j) = t2/betal;
		tt2(j) = 0.0;
		tt4    = (tt1(j)-1.0)*log(t2)-t2/betal-tt1(j)*log(betal);
		for (i=4;i<=30;i++)
		{
			t5        = double(i);
			tt3(j,i)  = (tt1(j)-1.0)*log(t5)-t5/betal-tt1(j)*log(betal)-tt4;
			if (tt3(j,i) < -1.0e30) //discontinuity here.
				tt3(j,i) = -1.0e30;
			tt3(j,i)  = mfexp(tt3(j,i));
			tt2(j)    = tt2(j) + tt3(j,i);
		}
	}
	
	for (j=1;j<=nbins;j++)
	{
		for (i=4;i<=30;i++) 
			tt3(j,i) = tt3(j,i)/tt2(j);
	}
	
	for (j=1;j<=nbins;j++)
	{
		for (i=j;i<=j+10;i++)
		{
			ii = (i-j)*5;
			if (ii == 0)
			{
				ttt(j,i)=tt3(j,1)+tt3(j,2)+tt3(j,3)*0.5;
			}
			else
			{
				ttt(j,i) = tt3(j,ii-2)*0.5+tt3(j,ii-1)+tt3(j,ii)+tt3(j,ii+1)+tt3(j,ii+2)+tt3(j,ii+3)*0.5;
			}
		}
		for (i=nbins+1;i<=nbins+10;i++) 
			ttt(j,nbins) = ttt(j,nbins) + ttt(j,i);
	}
	
	for (j=1;j<=nbins;j++)
	{
		for (i=1; i<=nbins;i++) 
			len_len(j,i) = ttt(j,i);
	}

FUNCTION totalselectivity
	int t, j;
	//SJDM  logistic with two periods, these periods should be in the data file.
	for(t = commstyr;t <= commendyr+1; t++)
	{
		for(j = 1; j <= nbins; j++)
		{
			if(t<= 2004)
				totalselect(t,j)=1./(1.+mfexp(-c2*(double(length_bins(j))-mfexp(log_d2))));
				
			if(t>2004)
				totalselect(t,j)=1./(1.+mfexp(-c3*(double(length_bins(j))-mfexp(log_d3))));
		}
		maxsel_totsel=max(totalselect(t));			//SJDM Another minor discontinuity.
		totalselect(t)=totalselect(t)/maxsel_totsel;
	}

FUNCTION retselectivity
	int t,j;
	//SJDM Logistic selectivity with two periods, breaking at 2004.
	//TODO use init_bounded_number_vector for selex parameters, can set up a number of periods in a control file.
	for(t = commstyr;t <= commendyr+1; t++)
	{
		for(j = 1; j <= nbins; j++)
		{
			if(t<= 2004)
				retselect(t,j)=1./(1.+mfexp(-aa2*(double(length_bins(j))-mfexp(log_bb2))));
			
			if(t>2004)
				retselect(t,j)=1./(1.+mfexp(-aa3*(double(length_bins(j))-mfexp(log_bb3))));
		}
		maxsel_retsel =max(retselect(t));		//Discontinuous?
		retselect(t) = retselect(t)/maxsel_retsel;
	}

FUNCTION gselectivity
	int j;
	for(j = 1; j <= nbins; j++)
	{
		gselect(j)=1./(1.+mfexp(-aa4*(double(length_bins(j))-mfexp(log_bb4))));
	}
	maxsel_gsel =max(gselect);
	gselect = gselect/maxsel_gsel;

FUNCTION gmortality
	//SJDM  groundfish bycatch mortality rates
	//I'm not sure what is going on here, with all these different periods
	int t;
	gfmort_tem.initialize();
	mean_gfmort.initialize();
	for(t = commstyr+10;t <= commendyr; t++)
	{
		gfmort_tem(t)=mfexp(log_mean_Fground+Fground_dev(t));
	}
	//SJDM calculating average groundfish fishing mort in 1995 to 1999, why?
	//Apparently using the average f prior to 1995, awkward way to do this.
	temp7=0.;
	for(t = commstyr+10;t <= commstyr+14; t++)
	{
		temp7 += gfmort_tem(t);
	}
	mean_gfmort = temp7/((commstyr+14)-(commstyr+10)+1);
	//
	gfmort_final.initialize();
	for(t = commstyr;t <= commendyr; t++)
	{
		if(t<1995)
			{gfmort_final(t)=mean_gfmort;}
		else
			{gfmort_final(t)=mfexp(log_mean_Fground+Fground_dev(t));}
	}

FUNCTION dvariable drnorm(const dvector& o, const dvar_vector& p, const dvar_vector& sig2)
	/*
		This function calculates the robust normal likelihood.
		Written by Martell, intended to replace the redundant code in the
		objective_function routine.
	*/
	double twopi = 6.283185;
	dvar_vector t1 = 0.5*log( twopi*sig2 );
	dvar_vector resid2 = elem_div( square(o-p),2.*sig2 );
	dvar_vector t6 = t1 - log(mfexp( (-1.)*resid2 )+0.01);
	return( sum(t6) );

FUNCTION objective_function
	int t,j;
	//
	//Robust liklihood approach for length composition fit
	//
	like_retlencomp=0.;
	for(t = commstyr; t <= commendyr; t++)
	{
		like_retlencomp += like_wght(1)*drnorm( obsretlencomp(t),predretcatchlen(t),retsigmasquare(t) );
	/*for(j = 1; j <= nbins; j++)
		{
		//like_retlencomp  += like_wght(1)*(0.5*log(6.286*retsigmasquare(t,j))-(log(mfexp(-1.*square(obsretlencomp(t,j)-predretcatchlen(t,j))/(2.*retsigmasquare(t,j)))+0.01)));
		}
		cout<<tmp<<endl;
		exit(1);*/
	}
	//
	like_discdlencomp=0.;
	for(t = commstyr+5; t <= commendyr; t++)
	{
		like_discdlencomp += like_wght(2)*drnorm( obsdiscdlencomp(t),preddiscdcatchlen(t),discdsigmasquare(t) );
	/*for(j = 1; j <= nbins; j++)
		{
		like_discdlencomp += like_wght(2)*(0.5*log(6.286*discdsigmasquare(t,j))-(log(mfexp(-1.*square(obsdiscdlencomp(t,j)-preddiscdcatchlen(t,j))/(2.*discdsigmasquare(t,j)))+0.01)));
		}*/
	}
	//
	like_gdiscdlencomp=0.;
	for(t = commstyr+10; t <= commendyr; t++)
	{
		like_gdiscdlencomp += like_wght(3)*drnorm( obsgdiscdlencomp(t),predgdiscdcatchlen(t),gdiscdsigmasquare(t) );
	/*for(j = 1; j <= nbins; j++)
		{
		like_gdiscdlencomp += like_wght(3)*(0.5*log(6.286*gdiscdsigmasquare(t,j))-(log(mfexp(-1.*square(obsgdiscdlencomp(t,j)-predgdiscdcatchlen(t,j))/(2.*gdiscdsigmasquare(t,j)))+0.01)));
		}*/
	}
	//
	// Q?  Why do you have different weights on the retdcpue from 1998-2003 and 2004-2010?
	// Q?  Why do you add the small constant to eh retained cpue (0.001)?
	
	like_retcpue=0.;
	double tiny = 0.001;
	dvector obs = log(obstretdcpue(commstyr+13,commendyr-7)+tiny);
	dvar_vector pred = log(predtretcpue(commstyr+13,commendyr-7)+tiny);
	dvar_vector sig2 = retdcpuevar(commstyr+13,commendyr-7)+ret_var;
	like_retcpue = like_wght(4)*drnorm( obs,pred,sig2 );
	/*for(t = commstyr+13; t <= commendyr-7; t++)  // from 1998 to 2003
		{
		like_retcpue += like_wght(4)*(0.5*(log(6.286*(retdcpuevar(t)+ret_var))+square(log(obstretdcpue(t)+0.001)-log(predtretcpue(t)+0.001))/(retdcpuevar(t)+ret_var)));
		}*/
	//
	dvector obs2 = log(obstretdcpue(commstyr+19,commendyr)+tiny);
	dvar_vector pred2 = log(predtretcpue(commstyr+19,commendyr)+tiny);
	dvar_vector isig2 = retdcpuevar(commstyr+19,commendyr)+ret_var;
	like_retcpue += like_wght(18) * drnorm(obs2,pred2,isig2);
	/*for(t = commstyr+19; t <= commendyr; t++)  // from 2004 to 2010
		{
		like_retcpue += like_wght(18)*(0.5*(log(6.286*(retdcpuevar(t)+ret_var))+square(log(obstretdcpue(t)+0.001)-log(predtretcpue(t)+0.001))/(retdcpuevar(t)+ret_var)));
		}*/
	//
	like_discdcpue=0.;
	ivector iyr(commstyr+13,commendyr);
	iyr.fill_seqadd(1998,1);
	like_discdcpue = like_wght(5)*drnorm( log(obstdiscdcpue(iyr)+tiny),log(predtdiscdcpue(iyr)+tiny),discdcpuevar(iyr)+discd_var );
	/*for(t = commstyr+13; t <= commendyr; t++)   // start from 1998
		{
		like_discdcpue += like_wght(5)*(0.5*(log(6.286*(discdcpuevar(t)+discd_var))+square(log(obstdiscdcpue(t)+0.001)-log(predtdiscdcpue(t)+0.001))/(discdcpuevar(t)+discd_var)));
		}*/
	//
	like_retdcatchB=0.;
	like_retdcatchB = like_wght(6) * norm2( log(obsretdcatchB+tiny)-log(predretdcatchB+tiny) );
	/*for(t = commstyr; t <= commendyr; t++)
		{
		like_retdcatchB += like_wght(6)*square(log(obsretdcatchB(t)+0.001)-log(predretdcatchB(t)+0.001));   // weight
		}*/
	//
	like_discdcatchB=0.;
	ivector jyr(commstyr+13,commendyr);
	jyr.fill_seqadd(commstyr+13,1);
	like_discdcatchB = like_wght(7) * norm2( log(obsdiscdcatchB(jyr)+tiny)-log(preddiscdcatchB(jyr)+tiny) );
	/*for(t = commstyr+13; t <= commendyr; t++)
		{
		like_discdcatchB += like_wght(7)*square(log(obsdiscdcatchB(t)+0.001)-log(preddiscdcatchB(t)+0.001));   // weight
		}*/
	//
	like_gdiscdcatchB=0.;
	ivector kyr(commstyr+10,commendyr);
	kyr.fill_seqadd(commstyr+10,1);
	like_gdiscdcatchB = like_wght(8) * norm2( log(obsgdiscdcatchB(kyr)+tiny)-log(predgdiscdcatchB(kyr)+tiny) );
	/*for(t = commstyr+10; t <= commendyr; t++)
		{
		like_gdiscdcatchB += like_wght(8)*square(log(obsgdiscdcatchB(t)+0.001)-log(predgdiscdcatchB(t)+0.001));   // weight
		}*/
	//
	like_rec_dev=0.;
	if(active(rec_dev))
	{
	for(t = commstyr+1; t <= commendyr+1; t++)
	{
	like_rec_dev += like_wght(9)*square(rec_dev(t));   //
	}
	}
	//
	like_F=0.;
	if(active(Fpot_dev))
	{
	for(t = commstyr; t <= commendyr; t++)
	{
	like_F += like_wght(10)*square(Fpot_dev(t));  //
	}
	}
	//
	like_gF=0.;
	if(active(Fground_dev))
	{
	for(t = commstyr+10; t <= commendyr; t++)
	{
	like_gF += like_wght(11)*square(Fground_dev(t));  //
	}
	}
	//
	//     like_Mpenalty= 0.;
	//    like_Mpenalty= like_wght(12)*square(M-0.18);    //0.18
	//
	//    like_initialLF = 0.;                           // 3rd order smoothness penalty for likely dome shaped first size distribution
	//   for(j = 1; j <= nbins-3; j++)
	//     {
	//   like_initialLF += like_wght(13)*(square(log(first_propN(j+3))-3.*log(first_propN(j+2))+3.*log(first_propN(j+1))-log(first_propN(j)))
	//                                   +square(log(first_propO(j+3))-3.*log(first_propO(j+2))+3.*log(first_propO(j+1))-log(first_propO(j))));
	//     }
	//
	like_qq=0.;
	if(active(qq_dev))
	{
	for(t = 1; t <= 2; t++)
	{
	like_qq += like_wght(14)*square(qq_dev(t));  //
	}
	}
	//
	like_molt50=like_wght(15)*(square(mfexp(log_b)-139.));  //There is no oldshell newshell data to estimate a & log_b
	// From tagging data
	like_beta=like_wght(16)*(square(mfexp(log_beta)-0.578));
	// Ratio of q pre and post rationalization
	like_meanq=like_wght(17)*(square(log(mfexp(logq3)/mfexp(logq2))-log(1.929204537)));
	//
	if(last_phase())
	{like_finalF =1.*square(Ftemp-final_F);}
	//
	f = like_retlencomp+like_discdlencomp+like_gdiscdlencomp+like_retcpue+like_discdcpue+like_retdcatchB+like_discdcatchB
	+like_gdiscdcatchB+like_rec_dev+like_F+like_gF+like_qq+like_molt50+like_beta+like_meanq+like_finalF;



REPORT_SECTION
	REPORT(mat_biomass);

TOP_OF_MAIN_SECTION
	time(&start);
	arrmblsize = 50000000;
	gradient_structure::set_GRADSTACK_BUFFER_SIZE(1.e7);
	gradient_structure::set_CMPDIF_BUFFER_SIZE(1.e7);
	gradient_structure::set_MAX_NVAR_OFFSET(5000);
	gradient_structure::set_NUM_DEPENDENT_VARIABLES(5000);
 

GLOBALS_SECTION
	/**
	\def REPORT(object)
	Prints name and value of \a object on ADMB report %ofstream file.
	*/
	#undef REPORT
	#define REPORT(object) report << #object "\n" << object << endl;

	#include <admodel.h>
	#include <time.h>
	//#include <stats.cxx>
	//#include <martool.cxx>
	time_t start,finish;
	long hour,minute,second;
	double elapsed_time;
	
FINAL_SECTION
	time(&finish);
	elapsed_time=difftime(finish,start);
	hour=long(elapsed_time)/3600;
	minute=long(elapsed_time)%3600/60;
	second=(long(elapsed_time)%3600)%60;
	cout<<endl<<endl<<"*******************************************"<<endl;
	cout<<"--Start time: "<<ctime(&start)<<endl;
	cout<<"--Finish time: "<<ctime(&finish)<<endl;
	cout<<"--Runtime: ";
	cout<<hour<<" hours, "<<minute<<" minutes, "<<second<<" seconds"<<endl;
	cout<<"*******************************************"<<endl;

